Critical scaling and aging in cooling systems near the jamming transition 
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We conduct athermal simulations of freely-cooling, viscous soft spheres around the jamming tran- 
sition density and find evidence for a growing length ^(t) that governs relaxation to mechanical 
equilibrium, ^(t) is manifest in both the velocity correlation function, and the spatial correlations 
in a scalar measure of local force balance which we define. Data for different densities (f) can be 
collapsed onto two master curves by scaling ^(t) and t by powers of — (;/>j|, indicative of critical 
scaling. Furthermore, particle transport for (f) > exhibits aging and superdiffusion similar to a 
range of soft matter experiments, suggesting a common origin. Finally, we explain how ^{t) at late 
times maps onto known behavior away from (j>j. 
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Discontinuous phases such as granular media, foams 
and emulsions often exhibit a jamming transition from a 
fluid to an amorphous solid as the volume fraction is in- 
creased This transition is characterised by the onset 
of elastic response, deriving from the deformation of par- 
ticles or droplets within a system-spanning contact net- 
work that inevitably arises when excluded volume con- 
straints cannot be satisfied. Numerical studies of model 
athermal systems have revealed similarities to continuous 
phase transitions: Static quantities such as elastic mod- 
uli vanish algebraically as thejamming volume fraction 
is approached from above [3| , accompanied by diverg- 
ing length scales in the linear response 3]. Furthermore, 
non-linear rheology under continuous shear has revealed 
a diverging viscosity as (/) ^ preconfiguring a finite 
yield stress above 0, H, [1] . The fiow curves can be 
made to collapse onto two master curves, one for above 
the transition and one below, when the stress and strain 
rates are scaled by powers of |0 — 0j| [1,0], analogous to 
critical scaling functions 

However, these static or driven systems mask the relax- 
ation mechanisms by which packings reduce unbalanced 
forces, preventing comparison with the rich phenomenol- 
ogy of relaxation in non-equilibrium systems. The noted 
similarities with criticality suggest we look there first: 
Systems quenched to their critical temperature attain lo- 
cal equilibrium on wavelengths shorter than a correla- 
tion length that grows algebraically with time Q. 
This growing length is at the root of the aging of cor- 
relation and response functions, which include factors of 
the form .g[^(t„ + t)/^(tw)] for a time since quench 
and a lag time t. After suitable normalization the func- 
tions g{x) are universal, belonging to a small number of 
dynamic universality classes. Aging is also present in 
(non-critical) glassy systems, but ^(i) is either implicit 
or has a different physical interpretation such as domain 
size 0. It is not known if any of this phenomenology 
carries over to athermal systems at or near <j)j. 

Here we numerically investigate non-linear relaxation 



in freely-cooling soft sphere systems for a range of 4> 
spanning j . We find evidence for a growing length scale 
^(t) corresponding to the relaxation towards mechanical 
equilibrium, i.e. force balance. In analogy with critical 
scaling, ^(f) for different (j) can be collapsed onto (j) < 
and (j) > (j)j master curves by scaling ^ and t by suitable 
powers of \(f> — (j>j\. From this scaling we also infer an 
unjamming time that diverges as — > For (j) > (jjj, 
£^{t) ^ t for all t attributable to elastic propagation, and 
the particle transport properties obey aging similar to ex- 
periments on a range of soft-matter systems 13, U , l3] , 
suggesting a common origin. Finally, for late times we 
show how ^(t) maps onto known behavior away from (pj. 

Model. — We consider viscous soft discs in a square sim- 
ulation cell of dimensions L x L with periodic bound- 
ary conditions. To reduce ordering effects, the parti- 
cle diameters d are uniformly distributed over the range 
[0.7(d), 1.3((i)] with (d) the mean. Particles have equal 
mass density. Overlapping particles a and P with cen- 
ters separated by a distance < ^{d" + d^) in- 
teract in two ways: (i) A linear repulsive force f^^ = 
— 2R°'^/{d°' + d^)] acting along the line of centers, 
and (a) a, viscous damping f™ = 7y(v" — v'') which acts 
to reduce the relative velocity between a and /3. Parti- 
cles making no contacts move ballistically. Note that for 
simplicity there is no distinction between tangential and 
normal velocities in the damping term. 

A structureless initial configuration is constructed by 
depositing particles uniformly over the simulation cell to 
give the required area fraction (j) = L~^^^Tr{d" /2)'^. 
Each velocity is initially zero, so the initial energy reser- 
voir is provided by overlapping particles and varies as (p"^ . 
This variation with (j) is weak compared to the much more 
rapidly varying quantities discussed below, so no cru- 
cial dependence on the choice of initial conditions is ex- 
pected. Particles were iterated using the velocity Verlet 
algorithm which includes particle inertia [l3| . The simu- 
lations were continued until the pressure either changed 
by less than one part in 10^ over a time interval corre- 
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spending to 20% of the total, or dropped below some pre- 
defined value. For each (j), the system size L was system- 
atically increased until the pressure and potential energy 
agreed over the entire time range, to within error bars. 
We observed that larger system sizes were required as 
(j)j was approached. Below we normalize lengths by (d), 
times by to = \J {d) (m) / ^ and the damping coefficient r] 
by rjo = with (m) the mean particle mass. 




FIG. 1: The static structure factor S{q,t) versus time for 
(a) (p = 0.7 and (b) (j) = 1.0. The time axis has been scaled 
by to and the q-axes by = 2n/{d). 

Results. — All results below are for a damping coeffi- 
cient r]/riQ = 0.04, which corresponds to a coefficient of 
restitution « 0.88 in the dilute limit when there are no 
multiple contacts. An overview of system evolution is 
provided in Fig. [H which shows the static structure fac- 
tor S{q,t) = N~^{p{c[,t)p{—q,t)) for high and low (f>, 
with N the number of particles and p(q, t) the spatially 
Fourier-transformed number density (here and through- 
out (• • • ) denotes averaging over different initial config- 
urations). Significant 0-dependence is seen only at late 
times, with the high density system displaying a weak 
signature of static large-length structure, in contrast to 
low densities where a peak emerges and grows in height 
and moves to smaller q with time. This corresponds to 



the cluster coarsening regime to be discussed later. 

For all (j) a local minimum-local maximum pairing 
emerges at short times and moves to lower q as the system 
evolves, approximately as ^ t~^ for small t. This struc- 
tural signature of a linearly growing length is also evident 
in a dynamic length extracted from the same-time veloc- 
ity correlations Cvv(r = |x"-x'^|,t) = Af{J2a,0^i^°' ^t) ■ 
v{x.^,t)), normalized so that Cvv{0,t) = 1. Following 
Olsson and Teitel [3| we identify the characteristic veloc- 
ity correlation length £,v{t) with the global minimum of 
Cvv • The growth of £,v{t) for different (j), and examples 
of Cm , are given in insets to Fig. [21 

The data for all (p can be collapsed onto two mas- 
ter curves by scaling £,v{t) by {(j) — and t by 
|(^ - with ly = 0.57 ± 0.05, e = 0.6 ± 0.05 and 
= 0.843 ±0.001, as demonstrated in Fig.H Note that 
we also scale t by a non-critical factor (f>^/^ to improve 
collapse at small times, in the spirit of corrections to scal- 
ing [7|, but this does not alter the exponents. The expo- 
nent V is consistent with that already found for steady 
flow but as for that protocol, we cannot achieve rea- 
sonable collapse using the exponent w 0.5 for the diverg- 
ing length in linear response 0], and conclude these two 
lengths are unrelated. Indeed, for all (j) > , ^v{t) t 
and hence diverges with time, whereas the linear response 
lengths only diverge for </> ^ (Convergence of e.g. 
pressure with system size is achieved as long as the dy- 
namics and statics have decoupled by the time ~ L/2). 
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FIG. 2: (Color online) Collapse of the velocity correlation 
lengths ^v(t) after scaling by \4> — 4>j\~'' and t hy\<j) — (j)j\~'^, 
with u = 0.57, e = 0.6 and 0j = 0.843. The upper line (blue, 
solid symbols) corresponds to </> > (jij, the lower line (red, 
open symbols) to (j) < (f)j. The rightmost points of data sets 
close to (j>j have been indicated. Dashed lines correspond to 
(X t and oc 1 — e~*^*^ and are intended to guide the 
eye. (Inset, lower right) PrecoUapsed data, left to right for 
decreasing <j). Each data set is truncated when ^v(t) ~ L/2. 
(Inset, upper left) Example of Cw for (j) = 0.92. 

Instead we propose that {t) corresponds to the length 
over which the system approaches mechanical equilib- 
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rium, i.e. force balance. To test this hypothesis, we 
assign to each particle a the scalar quantity = 1 — 
I S/3 f"''!/ X]/3 |f"^l> where the sums are over all particles 
(3 in contact with a and ("^ is the elastic-only component 
of the corresponding interaction force. Note that f "'^ 
is the resultant elastic force and |f"'^| a normalization 
factor, so higher ip means more balanced forces, with per- 
fect balance aX = 1. Proceeding as before, we measure 
spatial correlations C^^{r,t) in ?/'(x") = ip" and extract 
a characteristic length (t) , as described in Fig. [3l As 
with the velocity correlation length, initially grows lin- 
early and obeys the same scaling with |0 — as for ^v, 
confirming they both reflect the same relaxation process, 
although is statistically noisy and the reliable time 
window correspondingly smaller. Nonetheless this con- 
firms that mechanical equilibrium is reached on a growing 
length -^vW -^e^W- 
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FIG. 3: (Color online) Characteristic length for force bal- 
ance £,^{t) versus time under the same scaling as Fig. [2] with 
the same notation. The dashed line is £,^{t) oc t. (In- 
set) C^^{r,t) = N{ip{Q)ip{r) — at different times for 
4> = 0.92. ^^{t) is identified with the rightmost maximum. 

The master curves themselves give insights into the 
physics underlying relaxation. For short times £,v{t) oc t 
for all (j>, which, since there is no ballistic transport of 
particles (see below)., must be due to elastic propagation 
through multi-particle contact networks. Inverting the 
axes scahng gives a speed of sound c oc — which 
is finite at when e = v, which is within error bars. 
Below (j)j, £,v{t) also initially grows linearly before satu- 
rating at a finite value, at which point the data collapse 
fails and cluster coarsening (which is not controlled by 
4>j) begins. The allows us to define an 'unjamming' time 
when the plateau is reached, which by reversing the axes 
scaling is seen to diverge as ~ — . 

Aging. — If ^(i) controls the relaxation as claimed, the 
mean-squared displacement (MSD) /S.r'^{t^ + t,t^) = 
Ea |x"(tw + t) - x"(tw)|^ should be a function of 
^(^w + t)/^{tvj) [9]. For (j) > (j)j, ^{t) (X t so we expect 



Ar^(tw + t,t^) ^ git/t^) corresponding to full aging. 
Conversely, ^{t) approaches a c/i-dependent constant for 
(j) < (j)j and critical aging should cease (although some 
other form of aging may recur deep in the cluster coars- 
ening regime). Examples of Ar^(iw + iw) above and 
below (j>j are given in Fig. 21 For cf) > (f>j the MSD takes 
a fixed form which systematically scales to lower ampli- 
tudes and later times with increasing t^,. By contrast, 
for (f) < (j>j no such systematic scaling is apparent over 
the available time window. 

Aging has been experimentally observed in a range of 
relaxing soft-matter systems 10|, lul, Il2| , and we specu- 



late that the underlying mechanism may be the same as 
in this athermal system. To test this we must first quan- 
tify the observed tw-scaling for (p > We first smooth 
the data by fitting each iw~curve to the 3-parameter fit 
M(tw)/{f + [t/T(iw )]""}, as shown in Fig.lHb). M(t„) 
describes the variation in the overall amplitude of par- 
ticle transport with t„, T(tw) is a relaxation time, and 
a is the early-time growth exponent. For large t^ both 
M(tw) and r(tw) are expected to scale algebraically with 
tvf, so we fit each to the form ^(1 -I- tv^/B)^^--''^ resp., to 
extract 6m and br- The exponents a, 6m and br for each 
(j) are plotted in Fig. [5l^a) . Since there is little apparent 
variation with 4>, we can improve the statistics by averag- 
ing over aU (P, giving a = 1.51 ± O.OI, 6m = -0.98 ± 0.02 
and br — 0.84 ± 0.05. This latter value would appear to 
suggest subaging rather than the full aging br = I ex- 
pected for a critical point Q , but we cannot yet rule out 
systematic errors due to the chosen fitting functions. 

It is now possible to compare these findings to the 
equivalent quantities measured in the experiments 
[nl, [13]. We find three areas of agreement: (i) Superdif- 
fusive particle transport Ar^ ~ t"" with a « 1.5, as in- 
ferred from the speckle decay in experiments [lO| and 
directly measured here [for t ^ t(<w)]; (H) Aging, with 
experimental br in the range 0.77 to 1.8, compared to 
br ~ 0.84 here, and (Hi) Convective decay of the scat- 
tering vector, which has been interpreted as evidence for 
the ballistic motion of elastic strain deformations through 
the material ll|, [lj| . Elastic waves are also present in 



our system, and while we do not refute this interpreta- 
tion of the data, it is interesting to note that we observe 
a linearly-growing length in this model, namely the cor- 
relation length ^(i) ~ t detailed above. We hypothesise 
that this may be the true origin of the experimentally- 
observed ballistic growth law. In this context, closer com- 
parison with experimental data would be desirable. 

Late times. — For (p > (pj, ^(t) becomes arbitrarily 
large and the system reaches global mechanical equilib- 
rium. The scaling properties of these static systems have 
been studied by non-inertial algorithms 0, 13] ■ For com- 
parison, in Table H] we give the corresponding quantities 
for states generated by our inertial algorithm. The agree- 
ment with previous results confirms the robustness of the 
static state to the preparation procedure. As our density 
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FIG. 4: (Color online) Mean squared displacement Ar (t^ + 
t, tw) for single runs at (a) <^ = 0.8 < <i>j and (b) = 0.88 > 
0J. In both cases, lines from top to bottom correspond to 
geometrically increasing in the range 0.4 < tw/to < 40. In 
(b) the smooth green lines are smoothing fits (see text). 
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FIG. 5: (Color online) ('a^ Transport and aging exponents ver- 
sus density. Smooth lines are the average over all (j). (b) Char- 
acteristic cluster size and velocity correlation length at 
late times for </> = 0.7. 



window is somewhat broad, extending to roughly 15% 
either side of c/ij, simple scaling cannot be assumed Q 
and a correction to scaling was included when fitting. 

For 4> < the system unjams and an unstable cou- 
pling between dissipation and density fluctuations leads 
to system-wide mass separation into clusters and voids. 
In this regime the small-q peak in S{q, t) is pronounced, 
allowing a characteristic cluster length £,c{t) to be ex- 
tracted, as plotted in Fig.[5]Jb). The velocity correlation 
length is also shown, and is similar to where their 
time windows overlap, demonstrating that £,{t) tracks 
cluster growth after unjamming. Over the available 
data window ^c{t) « ^o.26±o.03^ towards the lower end 
of quoted values for hard spheres [lB|, suggesting that 



TABLE I: Fits of dimensionless pressure, P.E., excess coor- 
dination number, relative particle overlap and shear modulus 
at t = oo to A{(j} — + B{(f) — (/ij)], with a correction to 

scaling oc B. Numbers in brackets denote uncertainty in the 
last digit. The bulk elastic modulus remains finite as <^ . 



Quantity 


Prefactor A 


0./ 


Exponent a 


Pressure 


0.9(1) 


0.8433(4) 


1.03(4) 


Potential Energy 


0.05(1) 


0.8434(4) 


1.98(9) 


z - z'"" 


4.3(4) 


0.8432(4) 


0.55(4) 


Particle overlap 


0.4(1) 


0.8433(5) 


1.03(4) 


G 1 Gaffinc 


4.8(4) 


0.844(1) 


0.47(3) 



subsequent evolution is described by granular gas theory. 

Conclusions. — The relaxation of freely cooling ather- 
mal systems considered here support the view that the 
model 'point-J' jamming transition is closer in nature 
to a continuous phase transition than a glass transition: 
While aging is common to both, the critical scaling of 
^(<) in Fig. El controlled by a single point is expected 
only for critical points. The correspondence is of course 
not complete, and as for continuous shear Q we expect 
the master curves to depend on the interaction potential. 
Nonetheless we believe a fundamental understanding of 
this important transition will best approached from the 
standpoint of critical point theory, and further modelling 
in this direction would be desirable. 
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